#ifndef PCPS_MATRIX_ALGEBRA_HEADER
#define PCPS_MATRIX_ALGEBRA_HEADER

#include <src/pcps.h>
#include <src/mantissa_exponent.h>

namespace pcps {

  //-------------------------------------------------------------------------------
  // pcps::SingularMatrixException -- An exception to be thrown when a matrix
  //                                  cannot be inverted
  //-------------------------------------------------------------------------------

  class SingularMatrixException : public std::runtime_error {

    public:

      // constructor from string
      SingularMatrixException(const std::string &str) : std::runtime_error(str) {}

  };

  //-------------------------------------------------------------------------------
  // pcps::matrix_inverse_lu -- Computes a matrix's inverse and determinant using
  //                            an LU decomposition.
  //-------------------------------------------------------------------------------

                                                 // exit status    array size      description
  template <class SCALAR>                        // -----------    ----------      ----------------------------------------------------------------
  void matrix_inverse_lu(const int n,            // unchanged                      dimension of the matrix
                         SCALAR & det,           // overwritten                    On exit, the matrix's determinant.
                         SCALAR * const A,       // overwritten      n*n           On entry, the matrix to be inverted.  On exit, the inverse.
                         SCALAR * const work,    // overwritten      n*n           work space for computing the inverse
                         int * const ipiv,       // overwritten      n             work space for pivoting during the LU decomposition
                         int * const inv_perm);  // overwritten      n             work space for inverting the LU pivoting

  template <class SCALAR> pcps::MantExp<SCALAR> matrix_determinant_lu(const int n, SCALAR * const A, int * const ipiv);

} // end namespace pcps

#endif
